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Abstract 

Significant progress in the development of efficient and fast algorithms for quantum 
chemical calculations has been made in the past two decades. The main focus has al- 
ways been the desire to be able to treat ever larger molecules or molecular assemblies — 
especially linear and sub-linear scaling techniques are devoted to the accomplishment 
of this goal. However, as many chemical reactions are rather local, they usually involve 
only a limited number of atoms so that models of about two hundred (or even less) 
atoms embedded in a suitable environment are sufficient to study their mechanisms. 
Thus, the system size does not need to be enlarged, but remains constant for reactions 
of this type that can be described by less than two hundred atoms. The question then 
arises how fast one can obtain the quantum chemical results. This question is not 
directly answered by linear-scaling techniques. In fact, ideas such as haptic quantum 
chemistry or interactive quantum chemistry require an immediate provision of quan- 
tum chemical information which demands the calculation of data in "real time". In 
this perspective, we aim at a definition of real-time quantum chemistry, explore its 
realm and eventually discuss applications in the field of haptic quantum chemistry. 
For the latter we elaborate whether a direct approach is possible by virtue of real-time 
quantum chemistry. 

Submitted for publication as a perspective article in Int. J. Quantum Chem. 
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1 Introduction 



The large number of possible nuclear configurations puts the curse of dimensionality on 
studies of the chemical reactivity of large molecular systems like metalloenzymes or transi- 
tion metal complexes in homogeneous catalysis. To find exactly those configurations which 
correspond to the reaction path searched for is very difficult. Although there exist several 
methods to sample the configuration space of the nuclear positions efficiently and in an unbi- 
ased way (see, e.g., Refs. [2j[3M|5]), a full ab initio treatment of the molecular system renders 
also these approaches very time demanding or even unfeasible. An ab initio treatment is, 
however, mandatory for the study of chemical reactivity involving the forming and breaking 
of chemical bonds. 

The overwhelming amount of data generated by computational methods calls for new ap- 
proaches to access it. In conventional methods, a change of nuclear coordinates {Ri} of a 
reactive or functional molecular system changes the potential energy^ 
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i.e., the electronic energy defined in the Born-Oppenheimer (BO) approximation. Once this 
energy has been calculated by a computer program, the results are collected and can be eval- 
uated. However, if the results were accessible immediately after applying the configurational 
change, the perception of the results of the calculations would be much more efficient since 
a real interaction between the scientist and the system under study could be established. 

In this perspective article we discuss a point of view on such a challenge that we may sum- 
marize under the term Real-time Quantum Chemistry. "Real-time" in this context means 
that the results of quantum chemical calculations are basically instantaneously available, 
i.e., on a time scale that a specific human sense would interpret as instantaneous. Quantum 
chemical information available in real time then offers a very efficient way to interactively 
study the reactivity of molecular systems. 

The structure of this work is as follows. We start with a definition of Real-time Quantum 
Chemistry and discuss its principles. Then, we review and evaluate existing techniques 
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to accelerate quantum chemical calculations with respect to their potential usage in Real- 
time Quantum Chemistry. Finally, we discuss the feasibility of a direct approach to Haptic 
Quantum ChemistrjHD which implements the ideas of Real-time Quantum Chemistry. In 
order to illustrate what can already be done with currently available quantum chemistry 
software we present some exploratory calculations. 

2 Principles of Real-time Quantum Chemistry 

A theoretical study of chemical reactivity requires the exploration of the potential energy 
surface. Hence, first the wave function of the molecular system has to be calculated, from 
which then the energy and the gradients can be obtained. Other properties like molecular 
orbital representations, the electron density, polarizabilities or partial charges could also be 
considered, but are secondary targets compared to energies and gradients. We thus define 

Definition Real-time Quantum Chemistry shall denote the very fast calcu- 
lation of the quantum mechanical response of a reactive molecular system in 
terms of the wave function and the corresponding energy and gradient due to 
a user-driven manipulation of the system's molecular structure such that the 
operator experiences the information flow without any time delay. 

Since finding the optimized wave function and calculating the energy and the gradient is 
essential to all reactivity studies, we define their instantaneous calculation as the core ob- 
jective of Real-time Quantum Chemistry. Obtaining additional properties in real time could 
then be considered an extension of (core) Real-time Quantum Chemistry. 

It is clear that the fast calculation of the core quantities is central also to ab initio molecular 
dynamics (AIMD) P Accordingly, the challenges which one faces in the context of Real-time 
Quantum Chemistry can also be found in the area of AIMD, where the fast calculation of 
the nuclear gradients (called ionic forces) is of paramount importance. Therefore, many 
techniques developed in AIMD to speed up the calculation of gradients are also important in 
a Real-time Quantum Chemistry framework. Although the forces are required as quickly as 
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possible, AIMD simulation results are interpreted after the simulation has been carried out. 
In Real-time Quantum Chemistry, however, we want to seamlessly merge the calculation and 
the perception of the calculated information. 

The seamless integration of complex quantum chemical information allows the interactive 
exploration of chemical reactivity. An approach which utilizes a force-feedback device as an 
input and an output device to transmit the quantum mechanical information to the operator 
is Haptic Quantum Chemistry (HQC) introduced by us in 200900 A force-feedback device as 
in HQC is also utilized in the interactive molecular dynamics framework,^ which implements 
a classical treatment of the forces in molecular systems. Another example is the recently 
presented semi-empirical interactive implementation for optimizing molecular structures of 
hydrocarbons during editing the structures'^] 



As it has been noted in Ref. [TT] a real-time visual experience requires at least ten frames 
per second (i.e. 100ms to update the system's state), if a frame represents a step in the 
shift of nuclear positions during structure optimization or reaction. In Haptic Quantum 
Chemistry, however, a smooth tactile experience requires an update rate of 1000 frames per 
second. Accordingly, the central question for real-time quantum chemical approaches is the 
following: Is it — in principle — possible to perform sufficiently accurate ab initio quantum 
chemical calculations in such a short time, i.e., in the millisecond range for a reasonably 
sized molecular model of a reactive system? 

The issue of accelerating quantum chemical calculations has been treated thoroughly before — 
mostly in the area of linear and sub-linear scaling techniquesli3i3 However, it is important to 
note that there is a fundamental difference between the problem formulated in the question 
above and the problem of achieving linear scaling. In Real-time Quantum Chemistry, we do 
not ask for methods which allow us to treat larger and larger systems but rather how fast 
we can perform a calculation for a particular system of constant size and target accuracy. 
This implies that one also has to focus on the prefactors and onsets of quantum chemical 
methods and not only on their overall scaling behavior. In other words, the actual execution 
time for a given system is the prime target. 

The reason why we focus on reactive systems with a (large but) constant number of atoms 
is that we concentrate on rather local chemical events, which are ubiquitous in chemistry. 



For example, they can be mediated by transition-metal centers as in enzymatic reactions or 
in homogeneous catalysis. Reactions at transition metals are usually restricted to a limited 
spatial region containing the metal center(s), its (their) ligand environment, the incoming 
reactant, co-factors, other reactants that may lead to side reactions and some proper model 
of the close environment in such a way that models on the order of 200 atoms are sufficient 
for a meaningful descriptionEMIS 

3 Methods to Accelerate Electronic Structure Calcu- 
lations 

Almost all acceleration techniques for electronic structure calculations developed so far aim 
at a linear scaling behavior especially for very large molecules (> 1000 atoms). A multitude 
of methods has been developed as discussed in many excellent reviews JHHHEH Here, we need 
to evaluate the existing methods from a different point of view. The overall goal is to decide, 
whether there is a certain lower limit of computation time for mid-sized molecular systems 
(50 — 200 atoms), and to identify approaches that are likely to be important for a real-time 
calculation of gradients and energies. 

In view of the fact that single-determinant models like Hartree-Fock theory and, most impor- 
tantly, density functional theory (DFT) are likely to be the best candidates for quantitative 
real-time reactivity exploration, we first need to consider the solution of the Roothaan-Hall 
equations. The two most important steps in a self-consistent solution of the Roothaan-Hall 
equations are the construction of the Fock matrix and the subsequent calculation of the 
density matrix. Of course the size of the basis set chosen is also very important, since it 
determines the size of the system. 

We skip the derivation of the one-particle mean-field Hartree-Fock theory^ and of Kohn- 
Sham (KS) density functional theorylilElH which lead both to an effective one-particle 
operator equation. For the sake of simplicity, the following equations are given for spin- 
restricted calculations. The generalization to the spin-unrestricted case is straightforward. 
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The Fock operator / for Hartree-Fock and Kohn-Sham models can be written ad^BGH 

N/2 

/(n) = h( ri ) + [ 2 tin) - 7 £«(ri)] + A v xc (n) , (2) 

a 

where 7 and A are two parameters that define the electronic structure model. For a Hartree- 
Fock calculation, one chooses A = and 7 = 1, whereas for a 'pure' Kohn-Sham calculation 
7 = and A = 1 with a non-vanishing v xc hold. With A = 1 and 7 6 [0, 1] hybrid approaches 
are described. 

The introduction of a finite basis set {4> u }, 



A 



Vi(r)=£)CU&,(r) (3) 

is the most convenient way to solve the spatial integro-differential self-consistent-field (SCF) 
equations that result when / is operating on an orbital ipi. This yields the well-known 
Roothaan-Hall equations, 

FC = SCe, (4) 



where C is the matrix of the expansion coefficients defined in Eq. (|3j), S is the overlap matrix 
and e is the matrix of the orbital energies q of the orbitals ipi. The elements of the Fock 
matrix F in the chosen basis {4> u } are given by 

F„ v = H c ;j c + J, v + 7 + A V™ . (5) 

Here, H core = T + V eN is the matrix representation of the one-electron operator, J is the 
two-electron Coulomb matrix operator, K is the two-electron exchange matrix operator and 
V xc is the two-electron exchange-correlation matrix operator. For the calculation of the 
matrices J and K the two-electron repulsion integrals are contracted with the elements of 
the density matrix P which is for real orbitals in closed-shell systems defined as 

N/2 

P^iv = 2 ^ C^ a C ua , (6) 

a 

and thus the reason for the self-consistent iterative solution of the Roothaan-Hall equation. 
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The elements of the matrix V xc are derived from the exchange-correlation density functional 
J 7 which can be a functional of the electron density p in the local density approximation, 
of the density and the gradient of the density Vp in the generalized gradient approximation 
(GGA), or of the density, the gradient of the density, and the kinetic energy density r 
in meta-GGA functionals. Since the matrix elements cannot be evaluated analytically the 
integration is usually performed on a mesh of grid points constructed from merged atomic 
grids 



M Nf M 



\v xc \ 
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where Nf 1 is the number of points {r a } in the grid of atom I and {w a } are the corresponding 
weights. Vi is the function that splits the molecular grid into atomic sub-grids. The electron 
density p at a grid-point r a is given by 



Eq. Q has to be solved iteratively, since F depends on the elements of C. After reaching self 
consistency the total energy of the system is calculated from the converged density matrix 
elements and the Fock matrix elements, 

£ BO = £ p„,>Hzr + \J2 P ^ t + ^] + E *w + ^ NN , (9) 

where ^/ NN represents the Coulombic pair interaction of all atomic nuclei. From this equation 
also the gradient V/ with respect to the nuclear coordinates {Ri} can be derived. Within 
a finite basis of, e.g., Cartesian Gaussian functions one obtains after a few rearrangements 
the expressiorJ^HILLZI 



Vi{pv\a\) - 7^V/(/iA|az/) 



v^bo = J2 P ^ (V/flJT) + \ £ P^Pxa 

\U> fj,u\cr 

+ J2 P*> ( V/v; x ;) - J2 ( v '^) + v^ NN (io) 
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It includes the Pulay forces 221 that are due to the origin dependence of atom-centered basis 
functions. In case of plane-wave basis functions these forces vanish and the expression can 
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be simplified. For a more compact notation an energy weighted density matrix Q, 

N/2 

Qvpu = 2 ^ tgC^gCva , (11) 
a 

has been introduced. The derivatives of the exchange-correlation matrix elements dV^/dRj 
depend on the specific nature of the exchange-correlation potential v xc . General formula- 
tions of the energy gradient can be found in Refs. |23j|24l|25j|26] for the local (spin) density 



approximation and in Ref. 27 for gradient-corrected functionals. 



3.1 Basis Sets 

Very important for any acceleration technique is the choice of the basis set in which the 
Kohn-Sham or Hartree-Fock orbitals are expanded. The larger the basis set the more 
accurate the calculations are, but also the more time they consume. Reducing their size is, 
therefore, a very seductive means to reduce computation time. In addition certain basis sets 
speed-up certain parts of the Fock matrix calculation but may have disadvantages in other 
parts. 

The two most widely employed explicit basis sets are linear combinations of Gauss-type 
atomic orbitals and plane waves. Plane waves are especially suited for the calculation of 
periodic and homogeneous systems whereas the Gauss-type orbitals are usually employed in 
molecular systems. Density matrices in Gaussian basis sets are therefore mostly sparse, i.e. 
band diagonal and thus promote linear-scaling techniques. 

Plane-wave basis sets are completely delocalized in the direct space but are very localized in 
reciprocal space, which is why they are usually applied in solid state physics. In molecular 
systems, however, a huge number of plane waves would be needed to obtain the same accuracy 
as with localized Gauss-type orbitals. On the other hand, some of the integrals can be 
calculated very fast by Fast Fourier Transforms (FFT) within a plane- wave basis set. Also 
the calculation of the forces on the nuclei is computationally less expensive since the basis 
functions are not position dependent. To increase the accuracy and limit the number of basis 
functions in molecular systems, hybrid codes combine Gaussians and plane waves JHHHH2] 
These developments allow for the calculation of molecules with up to a million atomsP^ 
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Pseudopotentiak@I53 can be employed for both types of basis sets in order to reduce the 
number basis functions. How many electrons are treated implicitly by the pseudopotential 
can, of course, be varied and determines the accuracy achieved. In plane-wave calculations 
pseudopotentials are crucial to avoid the description of the nodal structure of orbitals so that 
the number of plane-waves can be limited to a reasonably small value. Pseudopotentials are 
also employed for heavy elements to account for relativistic effects.^ 

A less common alternative to the classes of basis sets described above are the so-called 
waveletspHl^iMH which aim at combining the best of both worlds. They are localized in both 
direct and reciprocal space and the integrals can be calculated with very fast methods similar 
to Fast Fourier Transforms. 

It is obvious, that for Real-time Quantum Chemistry the number of basis functions needs to 
be as small as possible. Since high accuracy does not need to be the main goal, comparatively 
small basis sets can be employed. For the same reason pseudopotentials are beneficial. Since 
plane waves and basis sets of plane waves mixed with Gaussians have successfully been 
applied in AIMD calculations, they are also useful for Real-time Quantum Chemistry. 

3.2 Fock Matrix Calculation 

The calculation of the Fock matrix can be divided into two different parts. The calculation of 
the elements of the one-electron and two-electron matrix operators. The calculation of the 
exchange-correlation matrix elements occurring in KS-DFT calculations will be discussed 
separately. In almost all ab initio electronic structure calculations the Fock matrix construc- 
tion is the most time consuming part, because of the evaluation of all integrals. Therefore, 
many sophisticated algorithms have been devised to speed up their evaluation. 

The calculation of the one-electron Hamiltonian matrix requires the evaluation of 0(K 2 ) 
matrix elements, where K denotes the number of basis functions. The matrix elements 
consist of two terms, the kinetic energy term and the electron-nuclei interaction part. If 
localized basis functions are employed, the number of integrals for the kinetic term will 
increase only linearly with system sizeP^ For the Coulomb interaction between the electrons 
and the nuclei, multipole expansions can be applied to reduce the number of integrals and 
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achieve linear scaling.^ Although the evaluation of these matrix elements has a quadratic 
scaling, their contribution to the overall execution time is for mid-sized molecular systems 
very small and hence no problem for Real-time Quantum Chemistry. However, to even 
speed up this part, program parallelization can be employed efficiently as the integrals can 
be evaluated independently from one another. 

The two-electron repulsion integrals (ERI) in the calculation of the Coulomb- and the ex- 
change matrix elements are in principle four-index quantities and thus their evaluation for- 
mally requires 0(K 4 ) operations. The first and most effective way to reduce their compu- 
tation time is to discard all elements whose contribution is below a certain threshold. The 
most widely employed technique for such an integral screening is based on Schwarz-inequality 
integral estimatesJ^3l22H Another approach to screen for negligible matrix elements are the 
multipole-based integral estimates.^ They consider in addition the 1/R decay behavior 
between two charge distributions. In non-direct SCF calculations the integrals are pre- 
calculated and integral screening can only be done at the level of the integrals. In direct 
SCF methods also the density matrix elements with which the integrals are contracted can 
be taken into account for screening. Therefore, also large integrals can be neglected, if they 
have only a small weight assigned by the density matrix elements. 

The calculation of the surviving integrals can then be accelerated. The first general ap- 
proach is to fit the densities occurring in the Coulomb and exchange integrals with auxiliary 
basis sets, thus reducing the four index ERIs to two index quantities. This approach to 
fit the densities accelerates the evaluation of the Coulomb matrix elements by an order of 
magnitude. To obtain the auxiliary basis sets two different methods are employed. One 
is to determine them by a fitting procedurePUEHHEU which is performed before the actual 
electronic structure calculation for each atom type and basis set and results in additional 
libraries of auxiliary basis sets. This density-fitting approach is also known as the resolution 
of the identity (RI)J3HHZ1E11E11 Another method is to employ a Cholesky decomposition (CD) 
algorithm to determine the auxiliary basis set for each calculation separately F^l l 51 ^^^ CD 
methods are computationally more demanding, but generate an auxiliary basis set, which is 
"the best" for a given basis and they do not depend on pre-fitted parameters. Therefore, the 
error introduced by the technique is controllable and no extra auxiliary basis-set libraries 
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are needed. 

The calculation of Coulomb matrix elements can also be made more efficient by employing 
hierarchical multipole methods like the fast multipole method^B^l or the quantum chem- 
ical tree code.®' Here, the problem of calculating the Coulomb interaction between many 
electrons is approximated by a truncated multipole expansion. A combination of a multipole 
expansion method and auxiliary basis sets can be applied to reduce the computation time 
even further.®! 

For the exchange matrix elements special methods exist to obtain a linear-scaling behavior. 
The RI technique, which is very efficient for the Coulomb matrix elements, can also be 
applied for the exchange matrix elements not with the same efficiency thought In addition 
several methods ave been developed specifically for the exchange matrix elements. Examples 
are the 0(iV)-Exchange method®' (where N denotes the number of basis functions and is 
identical to our K), the local K algorithm,®! the LinK method®' or the auxiliary density 
matrix methods for Hartree-Fock-type exchange.®' 

In KS-DFT calculations the Fock matrix F contains additional matrix elements from the 
exchange-correlation functional. The matrix elements are evaluated by numerical integration 
of the functional derivative on a grid. The computational cost depends, of course, on the 
size of the grid. The overall molecular grid is constructed from atomic grids, which are 
merged to obtain the molecular grid. A very common choice is the Becke atomic partitioning 
scheme.®' Since the calculation for each atomic grid can be done independently the overall 
scaling behavior can be made linear ]1U[66][67][68] rp^g eX change-correlation matrix elements 
are usually considered as requiring a negligible amount of computation time compared with 
the Coulomb and exchange matrix elements. However, after accelerating the two latter by, 
e.g., RI methods their evaluation becomes the most time consuming part in the Fock matrix 
construction. 

Considering the potential of the above-mentioned techniques to accelerate the Fock matrix 
construction, certainly the fast calculation of ERIs is of prime importance. Based on the 
current state of the art a rigorous screening of the integrals and density-fitting techniques 
are the methods of choice. Since the integrals can be evaluated independently they are prone 
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for a parallelized evaluation, which will be discussed later. The exchange-correlation matrix 
elements require efficient numerical techniques for an accelerated calculation. Straightfor- 
wardly, the coarseness of the grids employed can be increased to decrease the computational 
effort within pre-defined accuracy bounds. 



3.3 Density- matrix Construction 

In conventional SCF-type electronic structure calculations the density matrix is calculated 
by a full diagonalization of the Fock matrix. This approach has the advantage that the 
complete eigenspace can be obtained, i.e., including all virtual orbitals. A disadvantage is, 
however, that the scaling with respect to system size is cubic (and can be easily prohibitive 
in plane- wave calculations). Due to the small pre-factor of the diagonalization procedure for 
Gauss-type orbitals it mainly poses a problem for systems with several thousands of atoms. 

Almost all methods discussed here have been developed for problems where the construc- 
tion of the Fock matrix is very fast and therefore the construction of the density matrix 
becomes the bottleneck of the calculations. An excellent review over methods, which avoid 



the diagonalization in semi-empirical calculations can be found in Ref. 69 An assessment of 
density-matrix methods for self-consistent-field calculations by comparing purification and 
minimization methods has recently been published by Rudberg et alW^ 

The diagonalization of the Fock matrix does not directly yield the density matrix but the 
coefficient matrix C which is then contracted to obtain the density matrix P. Although 
the density matrix is a local quantity due to the nearsightedness of the electrons,^ the 
coefficient matrix is not. Therefore, a great reduction of computation time would be achieved 
calculating the density matrix directly from of the Fock matrix. 

The so-called energy minimization techniques exploit the fact that the correct density matrix 
minimizes the expression tr[PF]. This minimization, however, has to be done under the 
constraints that the density matrix fulfills the idempotency condition and the trace condition, 
which requires that the trace of the density matrix yields the number of electrons. By 
contrast, the diagonalization procedures explicitly fulfill the idempotency condition. Li, 
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Nunes, and Vanderbilt proposed the functional 

E (P) = tr [(3P 2 - 2P 3 ) (F - fil)] , (12) 

as a method to compute the density matrix by implicitly fulfilling the idempotency condi- 
tion.^ Here, /x denotes the chemical potential and / the identity matrix. Although this 
functional requires the density matrix to be set up in an orthogonal basis, there is no major 
problem to obtain a formulation for non-orthogonal basis sets, for which the energy func- 
tional is modified to include also the overlap matrix.--- Independently, Daw proposed a similar 
approach in 1993 73 that applied additionally steepest-decent iterations to minimize the func- 
tional. Alternatively, a conjugate-gradient approach can be chosen, which is then called the 
conjugate-gradient density-matrix search method.^3^ Then, the Fock and the density ma- 
trix are transformed into an orthonormal basis. Another advantage of this approach is that 
the chemical potential does not need to be known in advance. Other techniques employ for 
example curvy steps^or include second derivatives^ to obtain faster convergence. All these 
methods employ a purification of the matrix proposed by McWeeny^ in order to fulfill the 
idempotency condition of the density matrix. The so-called sign matrix method8^M3 follow 
a different way to achieve this by expressing the density matrix in terms of the sign matrix 
function, which can be computed by iterative schemes. Also possible is the introduction of 
a penalty functional for violating the idempotency condition. 70 However, the penalty func- 
tional is difficult to employ since it cannot be minimized by standard methods. Linear scaling 
of these methods is obtained by weakening the idempotency condition through restricting 
the minimization to localized density matrices with density matrix elements corresponding 
to a distance larger than a certain threshold forced to be zero. Other methods, which are not 
that common in ab initio electronic structure calculations, shall just be mentioned here for 
the sake of completeness. Examples are the Fermi Operator ExpansionSHD and the Fermi 
Operator Projection^ method. 

Another possibility to minimize the energy without explicitly imposing the orthonormality 
condition, as it is done by full diagonalization, is to minimize with respect to orbitals only 
with an implicit orthonormalization constraint. The Orbital Minimization^ll^Hl^lll^ll or the 
Optimal Basis Density Matrix Minimization^^ methods are typical examples. However, 
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they are mainly applied in large tight-binding or semi-empirical calculations.^ A linear- 
scaling behavior can be achieved in these methods by carrying out the minimization with 
respect to localized orbitals^ (for example, by searching only over functions which non-zero 
values inside a specified region). This approach would not introduce any approximation if the 
localized orbitals could be obtained by a unitary transformation of the occupied eigenstates 
of the Roothaan-Hall equation. 

3.4 Acceleration of SCF Convergence 

Assuming that the build-up of the Fock matrix and the subsequent calculation of the density 
matrix are fast enough for Real-time Quantum Chemistry, the SCF procedure poses severe 
obstacles. The sheer number of SCF iterations strongly depends on the atomic configuration 
and is hardly predictable. Therefore, it cannot be guaranteed that the wave function and 
the gradients are available in time. From this point of view, for Real-time Quantum Chem- 
istry it would be desirable to have a method which completely avoids any iterative methods. 
However, in all true ab initio electronic structure calculation an iterative procedure is un- 
avoidable, since the Fock matrix depends on the density matrix elements. Consequently, 
to be able to obtain a real-time calculation, one has to focus on reducing the number of 
iterations to a minimum. Performing the structural manipulations in small steps is therefore 
the key for a working Real-time Quantum Chemistry implementation. 

The steps discussed in the previous sections (Fock-matrix assembly and density-matrix cal- 
culation) are both part of a single step in the self-consistent-field procedure. The convergence 
of SCF iterations strongly depends on the first guess for the molecular orbitals and on the 
nuclear configuration. In a reactivity study, however, it may happen that configurations of 
the atoms occur, for which the SCF procedure converges only slowly or not at all. It is thus 
of utmost importance to have methods at hand, which yield stable and fast converging SCF 
iterations. 

Direct inversion in the iterative subspace (DIIS) 9 3E1 is a widely employed technique to 
accelerate and stabilize SCF iterations. Error vectors from previous iterations are calculated 
and minimized in a least-squares sense. Accordingly, previous iterations are utilized to 
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extrapolate the Fock matrix in the next iteration. Quite closely related to the DIIS method 
are techniques called Fock matrix dynamics^ or electron density extrapolation^ which are 
common in the field of BO molecular dynamics. Instead of accelerating the SCF algorithm 
itself the whole single-point SCF calculation is accelerated by extrapolating the information 
from previous time steps of the simulation. It is assumed that in between two steps the 
nuclear coordinates change only little so that the molecular orbitals, Fock matrix elements, 
or the density do not differe much from one step to the next and hence are a good starting 
point for the electronic structure optimization for the new nuclear configuration. Although 
SCF acceleration schemes have a long history, significant improvements can still be made as 
highlighted by the augmented Roothaan-Hall method.^ 

The pseudo-diagonalization technique^ is based on the observation that only for the first 
and the last SCF iteration a full diagonalization is necessary. In between it is sufficient to 
eliminate all Fock matrix elements connecting the occupied and virtual molecular orbitals 
by unitary transformations. As a consequence, the diagonalization of the Fock matrix blocks 
corresponding to the virtual- virtual and the occupied-occupied molecular orbitals can be 
avoided, which greatly reduces the time of each SCF iteration. 

In cases of a too narrow gap between the highest occupied and the lowest unoccupied molec- 
ular orbital, slow convergence or even divergent SCF iterations can occur. In such cases 
level shifting^ can be applied to enlarge the gap and therefore avoid a mixing. With this 
procedure the problem of divergence, slow convergence, or oscillating behavior can be cured 
in many cases. 

In all techniques discussed above the construction of the density matrix and the optimization 
of the molecular orbitals in the SCF iterations were independent processes. However, for 
methods based on minimization of an energy functional for constructing the density ma- 
trix, methods have been developed that combine the density matrix optimization and the 
self-consistent-field iterations in one single optimization loop JHHHEUllooi The idea was de- 
veloped for the coupled electron-nuclei problem (Car-Parrinello molecular dynamics 97 ) and 
is therefore often called the "molecular dynamics" method, but it can also be applied in 
situations were the nuclei are kept fixed. The difference to conventional matrix diagonaliza- 
tion procedures to solve for the eigenstates is that the variational principle is applied in a 
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dynamic fashion and all eigenstates are determined simultaneously. The dynamic variables 
are here the coefficients of the basis functions with a fictitious electon mass. To exploit this 
methodology for a Real-time Quantum Chemistry implementation one could require that 
the manipulations are done in a continuous fashion, i.e., along a 'trajectory' with sufficiently 
small configuration-change steps. 

3.5 Gradient Calculation 

As outlined in the beginning of this section the calculation of the energy gradient with 
respect to the position of the nuclear coordinates involves contributions from each term in the 
electronic energy. An efficient force evaluation for large molecular systems in the framework 
of pure DFT has recently been proposed by Reine et alW^ by combining screening with 
a fast multipole method. In addition the calculation of the Coulomb contributions were 
accelerated by employing a density-fitting scheme^323 with auxiliary basis sets. There are 
also efficient gradient implementations available that do not not employ density fittingJ23i22l 

3.6 Subsystem Approaches 

To divide a molecular system under study into smaller subsystems is a key to the molecular- 
model approach discussed in the first two sections of this work and offers the possibility 
to reduce the computational effort further. As through a magnifying glass, the reactive 
part /region of a molecular assembly can be embedded in a spectator background and this 
magnifying lens can even be moved around in the whole systemP^H So-called combined 
quantum mechanics/molecular mechanics (QM/MM)G53E1!D approaches are widely employed 
to enable calculations of large enzymatic systems. In QM/MM the reactive part is treated 
quantum mechanically and the surrounding environment is modeled by classical force fields. 
By contrast, QM/QM methods apply the laws of quantum mechanics to all subsystem but 
may treat them with different methodologies!^ 

Depending on how the subsystems are embedded into each other a different level of accelera- 
tion can be achieved. Completely independent subsystem methods — also called Divide-and- 
Conquer (D&C) approaches — facilitate a massively parallelized calculation. In the original 
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formulation of the D&C approach the independently calculated density matrices of the sub- 
systems were merged to yield the full density matrixP^S! Although the computation of the 
density matrices is done independently, the surrounding of a fragment is accounted for by 
buffer regions around the fragment. This approach does not only accelerate the calculation 
of the density matrix but also the calculation of the Fock matrices in the SCF calcula- 
tions, which is the reason why this methodology is treated in this separate section and not 
together with other density matrix construction schemes. One difficulty for Real-time Quan- 
tum Chemistry is, however, that the calculation of the D&C force for the full system is not 
well defined. But there are cures available.'^' The partitioning can also be carried out for 
other quantities than the density-matrix. For instance, in the fragment molecular orbital 
theory the fragmentation is done at the level of the molecular orbitalsP^ For a recent and 



comprehensive review on fragmentation methods we refer to Ref. Ill 



Subsystem techniques allow for an in principle exact embedding of the subsystems thus re- 
taining only the approximations introduced by choosing different electronic-structure meth- 
ods for the subsystems. In the framework of density functional theory such methods have 
been proposed and are widely employed, for instance, to account for solvent effectsJi^EiSElIIIllIIIllinZIIIII 

To assess the potential acceleration gained by any subsystem approach, one has to consider 
two aspects: (i) how small can the subsystems be chosen and (ii) how are they connected 
to each other, if at all. This affects not only the computation but also the accuracy. But 
since we are interested in local phenomena, the fragmentation is often already implied by the 
structure of the molecular system itself. Clearly, the smaller the subsystems can be chosen 
and the less one has to account for embedding effects, the better for Real-time Quantum 
Chemistry. 

From the Real-time Quantum Chemistry point of view, Divide-and-Conquer and density 
embedding approaches are appealing, since they allow the greatest reduction in computation 
time if the whole system needs to be treated quantum mechanically. If, however, large parts 
of the molecular system are not directly involved in a reaction, but rather serve as a dielectric 
environment, then QM/MM methods are most suitable. 
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3.7 Technical Aspects: Parallelization and Special Hardware 



Essentially all of the most popular quantum chemistry codes can be run in parallel. How- 
ever, these software approaches usually have a rather high overhead making them not the 
best candidates for Real-time Quantum Chemistry applications. But electronic structure 
calculations have not only been accelerated by developing increasingly efficient algorithms 
but also by employing latest and even specialized hardware. 

Although graphical processing units (GPUs) were originally designed for the fast rendering 
of three dimensional graphics, they can be employed for the acceleration of quantum chem- 
ical calculations. In the past years the development and application of special algorithms 
for quantum chemical calculations that exploit the computing power of GPUs has gained 
considerable attentionl^lHH]I J^l l 123 IH2£lH^I This trend is due to a significant effort undertaken 
to program quantum chemical algorithms specifically designed for GPUs but also due to new 
graphic cards produced with a focus on scientific calculations. For example, the widely used 
GAMESS US packagers G r the Terachem prograrrpSl are ^ efficiently exploit the ad- 
vantages of GPUs for electronic structure calculations. Compared to calculations performed 
only on the central processing unit (CPU) of a computer they are able to achieve considerable 
speed-ups in the order of one magnitude. 11221 Also in the field of semi-empirical calculations 
GPUs have attracted some attention)^ 

Massive parallelization is also possible on processor architectures other than GPUs. For 
Kohn-Sham DFT calculation, for instance, the application of processors from ClearSpeed 
Technology Ltd. has been reported to accelerate electronic structure and gradient calcula- 
tionsli22ILL22l Besides these developments for specific processor types, there are also some 
more general considerations about how to exploit the advantages of emerging new processor 
types available in the literature ji^Uli-EI j± comprehensive overview of special processors and 
their potential for electronic structure calculation can be found in Ref. 133| In the field 



of classical molecular dynamics simulations the so-called ANTON processor developed by 
Shaw et alW^ is a successful attempt to build such a special purpose processor. Also for 
electronic structure calculations special processors have been designed; namely ERIC, the 
ERI Calculation specific chip-multiprocessor^' or the molecular orbital calculation specific 



embedded high performance computing (EHPC) systempSSl 

Almost all special processor types mentioned here accelerate the electronic structure calcu- 
lations by parallelization of the ERI calculation. This is the reason for the success of GPUs 
but also of high performance clusters like, for instance, the IBM Blue Gene series JHZlllssj 
The calculations of ERIs is in almost all quantum chemical calculation the bottleneck be- 
cause of their sheer number. Specialized hardware which allows for a massive parallelization 
can directly accelerate the calculations and not only improves the scaling behavior. Thus, 
exploiting these techniques is certainly imperative for Real-time Quantum Chemistry. 

4 Direct Haptic Quantum Chemistry 

After having elaborated on the available and future quantum chemical methods for Real-time 
Quantum Chemistry implementations, we shall now discuss their benefits for Haptic Quan- 
tum ChemistrjHll and subsequently discuss their capabilities in an out-of-the-box application 
presented in the next section. 

The concept and implementation of Haptic Quantum Chemistry as presented in Refs. |6[[7] 
is to employ a force-feedback device as depicted in Fig. [T] as an input and an output tool 
allowing for an intuitive manipulation of molecular structures while feeling the gradients on 
the manipulated atoms as forces. Such approaches are referred to as haptic enabled inter- 
active molecular visualizations systems in the literature.'^' There are a few such methods 
already available, but they only employ classical force fields to calculate the forces rendered 
by the device 0°l l 141 H^ll^l Thus, they prohibit the study of chemical reactions, which would 
require the ability tp form and break chemical bonds. 

In our current set-up, the haptic device is a pen-like pointer which allows the user to manip- 
ulate objects in a virtual reality framework and, at the same time, to physically experience 
a force feedback (cf. Fig. [I]). Hence, one is able to feel the curvature of the potential energy 
surface of the manipulated nuclear coordinates in the whole system. The visual presentation 
of the structure, the gradients on the atoms not manipulated with the device, the orbitals, 
the electron density and other properties allow one to perceive complex information in a very 
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intuitive way compared with the sole visual presentation. Hence, the haptic quantum chemi- 
cal approach immerses the scientist more into the scientific problem. Even deeper immersion 
can of course be achieved by employing 3D displays or 3D glasses and other techniques from 
the field of virtual reality. 




Figure 1: The haptic device in Haptic Quantum Chemistry showing how a bromine molecule 
is moved towards an ethene molecule. The tip of the haptic pointer corresponds to the 
bromine atom next to the ethene. 

Since the human haptic sense is much more sensitive than the visual sense, haptic devices 
usually have an update rate of about 1 — 4kHz. By contrast, to create the illusion of a 
smooth movement for the visual sense only 25 Hz are necessary. To circumvent this problem, 
in Haptic Quantum ChemistryEE so far the forces are not calculated directly but are obtained 
from interpolating single-point gradients {g}. The force fj acting on an atom / is calculated 
from the interpolated gradient g by 

fi = ~9i, (13) 
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where the gradient of the energy is given by 

g I = V I E BQ . (14) 

Here, V/ denotes the three Cartesian derivatives with respect the nuclear coordinate of 
atom I. The computationally inexpensive interpolation facilitates the very fast calculation 
of the forces necessary for the high update rate of the haptic device. The main drawback is, 
however, that a coarse-grained gradient field of the molecular system needs to be calculated 
in advance (though it can be refined during haptic exploration as more data points can be 
calculated in the background). 

The pre-calculation of single-point gradients becomes, however, more and more demanding if 
parts of the molecular structure are allowed to relax. The molecular structure at one specific 
position in the configuration space of the mobile part is no longer unique and depends on the 
trajectory leading to it. As a consequence, one has to design algorithms that keep track of 
the history of the actual haptic exploration run. However, employing a Real-time Quantum 
Chemistry framework would allow us to circumvent this problem as the quantum chemical 
data is always immediately available so that no history (trajectory) needs to be stored. We 
may call this approach Direct Haptic Quantum Chemistry (D-HQC). 

At first glance, the high update rate (more than 1 kHz) of the haptic device requires an 
update rate of a millisecond if the gradients from the electronic structure calculation were 
directly rendered by the device. But as it has been shown in the field of interactive fluid 
dynamics simulations, the servo loop of the haptic device and the simulation loop can be 
separated and can operate with different update ratesJ^ Both loops are connected by a 
shared memory from where the servo loop constantly reads the current force written by the 
simulation loop as soon as the new force is available. To reduce the occurring artifacts the 
stepwise force update is smoothed by a force filter. With this technique the update rate 
of the molecular system can be lowered below 100Hz. Accordingly, D-HQC would require 
that within a few hundred millisecond a new gradient on the manipulated atoms has to be 
available and the other atom positions have to be relaxed (Fig. [2|. 

This almost instantaneous relaxation of the whole structure is, however, not always wanted. 
By contrast, it might be even desired to be able to alter the molecular structure faster than 
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Figure 2: Program structure of a Direct Haptic Quantum Chemistry implementation. 

the system relaxes in order to simulate non-adiabatic changes. Therefore, the relaxation 
does not need to be instantaneous as long as the force change due to the relaxation can be 
rendered smoothly. For almost adiabatic changes the structure alterations have to be slowed 
down or done in small steps so that the structure of the whole system can relax. Altering the 
structure only in very small increments also speeds up the electronic structure optimizations, 
since the molecular orbitals change only very little. The program could force the user to 
perform only small changes (slow movements) by applying re-stalling forces on the haptic 
device. 

D-HQC can be described as probing the potential energy surface in the configuration sub- 
space spanned by the manipulated atoms of the molecular system. For a more formal de- 
scription of the force calculation in D-HQC the set of nuclei I is partitioned into nuclei 
controlled by the haptic device I' and the remaining nuclei I". The force on a nucleus I' 
is then calculated as the negative spatial derivative of the total energy, which is minimized 
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with respect to the basis set and the nuclear coordinates of the remaining nuclei {Rt"}. 

f p = -Vr min E tot [{</>*} ; {R ri Rp>}} (15) 

{4> t },{Rj„} 

with the total energy functional E tot is given by 

Etot JW] = (M{Ri',Ri"}) Hf * ({R r ,R r ,}j) . (16) 



5 Examples 

Many of the above-mentioned developments of linear and sub-linear scaling methods are 
available in standard program packages, though the ultimate package for real time quantum 
chemistry has not been developed yet — mostly for the reason that each package has been 
designed to serve a certain purpose. For example, huge molecular models have already been 
studied with the CP2KHHH program that employs mixed Gaussian and plane wave basis 
sets in AIMD simulations. Here, we choose the very efficient DFT modules of the Tur- 
bomole package^MEMI combined with our D-HQC setup to demonstrate that such studies 
are in reach for molecular models of relevant size. The calculations exploit density fitting, 
effective core potentials and small basis sets. It is clear that the resulting accuracy then 
does not necessarily live up to the current standard. However, this is also not decisive as 
a reaction pathway recorded during a D-HQC exploration can always be relaxed on a more 
accurate potential energy surface. 
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5.1 D-HQC for a bromine molecule approaching an ethene 




Figure 3: Exploratory calculation for the reaction of a bromine molecule (in red) with ethene. 
The distance from the closest bromine atom to the ethene is written above this bromine atom. 
Below it, the average time per single point in the optimization and the overall time until 
structure convergence are given. 



As depicted in Fig. 5T in this exemplary study a bromine molecule was pushed onto an 
ethene molecule to probe its reactivity. In this setup the position of the bromine atom 
closest to the ethene was dragged towards the ethene molecule so that the relaxation at 
each step had the bromine to ethene distance as a constraint. In this way a trajectory in 
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the subspace spanned by the position of the bromine atom was generated. The distance 
was incrementally decreased by 0.5 Bohr. The resulting distances in Angstroms are given in 



Fig. 5.1, which shows three exemplary points from the trajectory. In addition, the average 
time per structure relaxation step, i.e., the time for updating the system's structure, as well 
as the overall time needed to converge the structure (second line) are depicted. 

Following the results of the discussion about basis sets we chose for the bromine atoms the 
Stuttgart ECP-28-MWP pseudo potential^ and the def2-SV(P) basis set, for the carbon 
atom the def2-SV(P) basis setP^Sl and for the hydrogen atoms the STO-3G HONDO basis 
se gmm in order to obtain a very fast calculation of single-point energies and gradients. 
The calculations were performed employing the BP86 exchange-correlation density func- 
tionaU^ll 152 H 153 lE^M^I on a coarse numerical grid. In addition, also the resolution-of-identity 
technique (RI)P was applied to accelerate the calculations. 



The execution times in Fig. 5^ show that the time needed to update the system is almost 
constant during the trajectory, but the structure optimization time increases when the re- 
actants get closer, which indicates that it needs more steps to converge. As it was outlined 
before, the update rate is the important quantity for a real-time experience. The execution 
times per update step can be shortened by sampling the trajectory in smaller steps, which 
means that the SCF procedure can converge faster since the wave function does not change 
too much from step to step. 

The computations here were performed by running the individual Turbomole modules 
sequentially. Note that there is still room for improving the efficiency in terms of passing the 
information from one calculation to the next and avoid read/write accesses to the hard disk 
as these processes have not yet been optimized for a D-HQC implementation in standard 
programs like TURBOMOLE. 



5.2 D-HQC for an Sn2 reaction of fluoride with chloromethane 

Another example with a more pronounced effect of structural relaxation of the remaining 
atom positions is expected to demonstrate how this influences the systems update rate: a 
Sn2 type reaction. In this example a fluoride ion approaches a chloromethane molecule and 
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replaces the chlorine atom. The trajectory was recorded by moving the fluoride anion in 
steps of 0.1 Bohr towards the C atom. At each step the structure of the remaining atoms 
was optimized. The intermediate electronic structure optimizations were carried out again 
with small basis sets (def2-SVP 148 for the C, F, and CI atoms and STO-3G HONDCGHHo] 
for the H atoms) in combination with effective core potentials (ECP-10-MWBl 133 for CI and 
ECP-2-SDF for F) and BP86/RI on a coarse numerical grid. To reduce the number of SCF 
cycles in each geometry optimization cycle the molecular orbitals from the preceding point 
of the trajectory were taken. 




Figure 4: D-HQC exploration for the Sn2 reaction of fluoride (green) and chloromethane. 
The chlorine atom is printed in blue. The distance from the attacking fluorine anion to the 
central carbon atom is written below the fluorine atoms. Above the atoms the average time 
per single point in the optimization (in ms) and the overall time until structure convergence 
(in s) is given. 



In Figure 5.2 three intermediate points of the trajectory are shown. The average time per 
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electronic structure calculation in the structure optimization is almost constant and on the 
order of several hundred milliseconds. The overall time for the structure optimizer to reach 
convergence is, however, not constant and increases significantly when the attacking fluoride 
approaches its position in the transition state. The detailed timings show again that not 
the electronic structure optimization but the increased number of cycles in the structure 
optimization give rise to the increased overall time. As already discussed in the previous 
section this is not a severe issue as only the system updates have to be fast, which is this 
case in this example. The user would have to slow down the movement of the fluoride 
when approaching the C atom in order to obtain a reasonable minimum energy path. For a 
non-adiabatic simulation the movement can be faster although the remaining atoms of the 
system are not able to relax in time. 

6 Conclusions 

The possibility and necessity of obtaining the result of quantum chemical calculations in real 
time is beneficial in many respects for studying the reactivity of chemical systems and may 
change the way how quantum chemistry is done in the future. Not only the ever increasing 
amount of information provided by calculations but also the inherent complexity of chemical 
problems calls for new approaches like (Direct) Haptic Quantum Chemistry that require a 
Real-time Quantum Chemistry framework. The overview of currently available techniques 
to accelerate calculations provided here clearly showed that Real-time Quantum Chemistry 
is in reach and will be possible for relevant system sizes in the near future. 

The evaluation of existing algorithms and technology for Real-time Quantum Chemistry also 
demonstrated, however, that a paradigm change is needed. Almost all techniques presented 
here were not specifically designed to allow quantum chemical calculations of energies and 
gradients in real time. The aim of their development was the overall scaling behavior to 
allow the treatment of ever larger molecules or molecular systems. For Real-time Quantum 
Chemistry the focus needs to be on reducing the actual execution time for a fixed system size 
to around 100 ms. Although most of the currently available program packages in quantum 
chemistry have not been developed to allow the ultra-fast calculation of molecular systems 
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consisting of 100 — 200 atoms, the greatest potential for achieving considerable speed-up 
towards real time lies most probably in the activation of specialized hardware. Already in 
reach are calculations on GPUs which show a promising potential, but also completely new 
specialized hardware is desirable for Real-time Quantum Chemistry. 

The overwhelming amount and the complexity of the data generated by current quantum 
chemical calculations already limits their fast and intuitive evaluation. Haptic Quantum 
Chemistry offers a new approach to tackle this problem. The instantaneous availability 
of the wave functions and the gradients offered by Real-time Quantum Chemistry allows 
an even more convenient way of studying chemical reactivity, as we have discussed for the 
Direct Haptic Quantum Chemistry variant. The exploitation of the human haptic sense to 
present scientific data more intuitively is only a first step. A deeper immersion by employing 
techniques already developed in the field of virtual reality would be the ultimate goal of any 
development in this direction. 
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